Leave-one-site out LICA
# Before ComBat
# all site (IC maps)
# Components were ordered by variance explaining.
all=readMat('./2023_11_14_lica_n140_ctareadti_B1000nn_rsfc.mat')
all_nocb = do.call(rbind,lapply(1:5, function(ss)all$Morig40[[7]][[ss]][[1]] %*% diag(as.vector(all$Morig40[[3]][[ss]][[1]]))))[,order(all$Morig40[[15]], decreasing=TRUE)]
colnames(all_nocb) = paste('noCombat.All.IC',stringr::str_pad(1:40,2, pad = "0"), sep='')
all_cb = do.call(rbind,lapply(1:5, function(ss)all$Morig40.combat[[7]][[ss]][[1]] %*% diag(as.vector(all$Morig40.combat[[3]][[ss]][[1]]))))[,order(all$Morig40.combat[[15]], decreasing=TRUE)]
colnames(all_cb) = paste('Combat,All.IC',stringr::str_pad(1:40,2, pad = "0"), sep='')
# leave one out site (IC maps)
# Components were ordered by variance explaining.
site=readMat('2023_11_14_lica_n140_ctareadti_B1000nn_rsfc_by_site.mat')
xmat_no_cb=do.call(cbind,lapply(1:5,function(j){
tmpmat = do.call(rbind,lapply(1:5, function(ss)site$loso[[j]][[1]][[7]][[ss]][[1]] %*% diag(as.vector(site$loso[[j]][[1]][[3]][[ss]][[1]]))))[,order(site$loso[[j]][[1]][[15]], decreasing=TRUE)]
colnames(tmpmat) = paste('Site',j,'.IC',stringr::str_pad(1:40,2, pad = "0"), sep='')
return(tmpmat)}
))
xmat_cb=do.call(cbind,lapply(1:5,function(j){
tmpmat = do.call(rbind,lapply(1:5, function(ss)site$loso.combat[[j]][[1]][[7]][[ss]][[1]] %*% diag(as.vector(site$loso.combat[[j]][[1]][[3]][[ss]][[1]]))))[,order(site$loso.combat[[j]][[1]][[15]], decreasing=TRUE)]
colnames(tmpmat) = paste('Site',j,'.IC',stringr::str_pad(1:40,2, pad = "0"), sep='')
return(tmpmat)}
))
# Find the best matching components to the all site
cor_no_cb = abs(cor(all_nocb, xmat_no_cb, method='pearson'))
cor_no_cb = cor_no_cb[, apply(cor_no_cb, 2, function(s)all(is.na(s)))==FALSE]
cor_no_cb = cor_no_cb[apply(cor_no_cb, 1, function(s)all(is.na(s)))==FALSE, ]
cor_no_cb2 = cor_no_cb[, apply(cor_no_cb, 2, function(s)all(abs(s<0.5)))==FALSE]
cor_no_cb2 = cor_no_cb2[apply(cor_no_cb2, 1, function(s)all(abs(s<0.5)))==FALSE,]
cor_cb = abs(cor(all_cb, xmat_cb, method='pearson'))
cor_cb = cor_cb[, apply(cor_cb, 2, function(s)all(is.na(s)))==FALSE]
cor_cb = cor_cb[apply(cor_cb, 1, function(s)all(is.na(s)))==FALSE, ]
cor_cb2 = cor_cb[, apply(cor_cb, 2, function(s)all(abs(s<0.5)))==FALSE]
cor_cb2 = cor_cb2[apply(cor_cb2, 1, function(s)all(abs(s<0.5)))==FALSE,]
#apply(cor_no_cb.tmp,2, which.max)
colist = c(rep("#FFFFE5",50),COL1('YlOrBr', 100))
drawcorplot <- function(mat, threshold=0.5){
mat[abs(mat)<threshold]<-0
mat = mat[apply(mat, 1, function(s)all(abs(s<threshold)))==FALSE,]
mat = mat[,apply(mat, 2, function(s)all(abs(s<threshold)))==FALSE]
corrplot(t(rerow(t(mat))),col = colist, col.lim=c(0,1),is.corr = FALSE)
}
Spatial pattern matching
no combat
drawcorplot(cor_no_cb2[,grep('Site1.', colnames(cor_no_cb2))], threshold=0.5)
## [1] "16 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."

drawcorplot(cor_no_cb2[,grep('Site2.', colnames(cor_no_cb2))])
## [1] "23 row has all zeros, delete this column."

drawcorplot(cor_no_cb2[,grep('Site3.', colnames(cor_no_cb2))])
## [1] "24 row has all zeros, delete this column."

drawcorplot(cor_no_cb2[,grep('Site4.', colnames(cor_no_cb2))])
## [1] "18 row has all zeros, delete this column."
## [1] "20 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."

drawcorplot(cor_no_cb2[,grep('Site5.', colnames(cor_no_cb2))])
## [1] "17 row has all zeros, delete this column."
## [1] "19 row has all zeros, delete this column."
## [1] "22 row has all zeros, delete this column."
## [1] "22 row has all zeros, delete this column."

combat
drawcorplot((cor_cb2[,grep('Site1.', colnames(cor_cb2))]))
## [1] "20 row has all zeros, delete this column."

drawcorplot((cor_cb2[,grep('Site2.', colnames(cor_cb2))]))
## [1] "17 row has all zeros, delete this column."
## [1] "20 row has all zeros, delete this column."
## [1] "30 row has all zeros, delete this column."

drawcorplot((cor_cb2[,grep('Site3.', colnames(cor_cb2))]))
## [1] "15 row has all zeros, delete this column."

drawcorplot((cor_cb2[,grep('Site4.', colnames(cor_cb2))]))
## [1] "8 row has all zeros, delete this column."
## [1] "21 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."

drawcorplot((cor_cb2[,grep('Site5.', colnames(cor_cb2))]))
## [1] "17 row has all zeros, delete this column."
## [1] "18 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."

Agreement computation
computemat <- function(mat,ncomp=20){
tmp=rerow(mat)[1:ncomp,1:ncomp]
re=data.frame(diag.element = mean(diag(tmp)), offdiag.element = 1-mean(tmp[lower.tri(tmp, diag=FALSE)]) )
return(re)
}
dat=do.call(rbind, lapply(c(5,10,15,20,25),
function(j){
data.frame(site = 1:5, ncomp=j,
combat = do.call(rbind,lapply(paste('Site',1:5,'.',sep=''),function(str)computemat(t(cor_cb[,grep(str, colnames(cor_cb))]),ncomp=j))),
no.com=do.call(rbind,lapply(paste('Site',1:5,'.',sep=''),function(str)computemat(t(cor_no_cb[,grep(str, colnames(cor_no_cb))]),ncomp=j))))
}))
fig1=dat %>%
melt(., id=c('site','ncomp')) %>%
mutate(method = factor(substr(variable,1,6), levels=c('combat','no.com'), labels=c('ComBat','no-ComBat'))) %>%
mutate(element = factor(substr(variable,8,11), levels=c('diag','offd'), labels=c('mean(Diagonal Elements)','1-mean(Off-Diagnonal Elements)'))) %>%
ggplot(., aes(x=site, y=value, color=method, group=method)) +
geom_point() + facet_grid(element~ncomp) +
geom_line() + ylim(c(0.5,1)) +
ylab('Similarity') +
theme_bw() +
theme(legend.position='bottom')+
labs(caption='* Higher values indicates higher similarity between matched components.\n * The similarity measures range from 0 to 1.')
fig1

png('OCDglobal_n140_SiteSimilarity_B1000.png', width=2400, height=1800, res=300)
fig1
dev.off()
## quartz_off_screen
## 2
fig1_1=dat %>%
select(-combat.offdiag.element,-no.com.offdiag.element) %>%
melt(., id=c('site','ncomp')) %>%
mutate(method = factor(substr(variable,1,6), levels=c('combat','no.com'), labels=c('ComBat','no-ComBat'))) %>%
ggplot(., aes(x=site, y=value, color=method, group=method)) +
geom_point() + facet_grid(~ncomp) +
geom_line() + ylim(c(0.5,1)) +
ylab('Similarity') +
theme_bw() +
theme(legend.position='bottom')+
labs(caption='* Higher values indicates higher similarity between matched components.\n * The similarity measures range from 0 to 1.')
fig1_1

png('OCDglobal_n140_SiteSimilarity_diag_B1000.png', width=2400, height=1200, res=300)
fig1_1
dev.off()
## quartz_off_screen
## 2
Combat vs. no combat comparison
combatmat=abs(cor(all_cb,all_nocb))
combatmat = combatmat[, apply(combatmat, 2, function(s)all(is.na(s))) == FALSE]
combatmat = combatmat[apply(combatmat, 1, function(s)all(is.na(s))) == FALSE, ]
thresh=0.5
combatmat2=combatmat;combatmat2[abs(combatmat)<thresh]<-0
combatmat2 = combatmat2[apply(combatmat2,1,function(ss)all(ss<thresh))==FALSE, apply(combatmat2,2,function(ss)all(ss< thresh))==FALSE]
#combatmat2 = combatmat2[apply(combatmat2,1,function(ss)all(ss<thresh))==FALSE, apply(combatmat2,2,function(ss)all(ss< thresh))==FALSE]
#combatmat2 = combatmat2[, apply(combatmat2,2,function(ss)all(ss<thresh))==FALSE]
#corrplot(rerow(combatmat2),col = colist, col.lim=c(0,1),is.corr = FALSE)
#corrplot(rerow(t(rerow(combatmat2))),col = COL1('YlOrBr', 200), col.lim=c(0,1),is.corr = FALSE)
#corrplot((rerow(combatmat2)),col = COL1('YlOrBr', 200), col.lim=c(0,1),is.corr = FALSE)
thresh=0.5
rerowed.mat = rerow(t(combatmat2))[c(1:31,33,32),]
## [1] "8 row has all zeros, delete this column."
corrplot(t(rerowed.mat),col = colist, col.lim=c(0,1),is.corr = FALSE)

pp=min(nrow(rerowed.mat),ncol(rerowed.mat));
matched.dat = data.frame(combatics = colnames(rerowed.mat)[1:pp],nocombatics = row.names(rerowed.mat)[1:pp],
matchingcor = diag(rerowed.mat[1:pp,1:pp])) %>%
filter(matchingcor>0.5)
knitr::kable(matched.dat) %>% kableExtra::kable_classic(full=F)
|
combatics
|
nocombatics
|
matchingcor
|
|
Combat,All.IC01
|
noCombat.All.IC01
|
0.9954084
|
|
Combat,All.IC02
|
noCombat.All.IC02
|
0.9987395
|
|
Combat,All.IC03
|
noCombat.All.IC03
|
0.9875149
|
|
Combat,All.IC04
|
noCombat.All.IC04
|
0.9967404
|
|
Combat,All.IC05
|
noCombat.All.IC06
|
0.9919961
|
|
Combat,All.IC06
|
noCombat.All.IC05
|
0.8671488
|
|
Combat,All.IC07
|
noCombat.All.IC09
|
0.9887035
|
|
Combat,All.IC09
|
noCombat.All.IC17
|
0.7780303
|
|
Combat,All.IC10
|
noCombat.All.IC12
|
0.9805605
|
|
Combat,All.IC11
|
noCombat.All.IC11
|
0.9772330
|
|
Combat,All.IC12
|
noCombat.All.IC25
|
0.9079190
|
|
Combat,All.IC13
|
noCombat.All.IC10
|
0.9648221
|
|
Combat,All.IC14
|
noCombat.All.IC13
|
0.9941005
|
|
Combat,All.IC15
|
noCombat.All.IC18
|
0.9505280
|
|
Combat,All.IC16
|
noCombat.All.IC20
|
0.9584342
|
|
Combat,All.IC17
|
noCombat.All.IC16
|
0.9743688
|
|
Combat,All.IC18
|
noCombat.All.IC22
|
0.9728010
|
|
Combat,All.IC19
|
noCombat.All.IC26
|
0.9731191
|
|
Combat,All.IC20
|
noCombat.All.IC30
|
0.6589200
|
|
Combat,All.IC21
|
noCombat.All.IC23
|
0.9628430
|
|
Combat,All.IC22
|
noCombat.All.IC24
|
0.9302998
|
|
Combat,All.IC23
|
noCombat.All.IC15
|
0.8817485
|
|
Combat,All.IC24
|
noCombat.All.IC14
|
0.9196102
|
|
Combat,All.IC25
|
noCombat.All.IC21
|
0.9062225
|
|
Combat,All.IC26
|
noCombat.All.IC29
|
0.8700443
|
|
Combat,All.IC27
|
noCombat.All.IC19
|
0.8775439
|
|
Combat,All.IC28
|
noCombat.All.IC28
|
0.9717423
|
|
Combat,All.IC29
|
noCombat.All.IC32
|
0.8331042
|
|
Combat,All.IC30
|
noCombat.All.IC33
|
0.9382460
|
|
Combat,All.IC31
|
noCombat.All.IC31
|
0.9338284
|
|
Combat,All.IC32
|
noCombat.All.IC35
|
0.9173853
|
|
Combat,All.IC33
|
noCombat.All.IC34
|
0.9226888
|
|
Combat,All.IC34
|
noCombat.All.IC08
|
0.5759699
|
#corrplot((rerow((combatmat2))),col = colist, col.lim=c(0,1),is.corr = FALSE)
#corrplot(rerow(t(rerow((combatmat2)))),col = colist, col.lim=c(0,1),is.corr = FALSE)
#corrplot(rerow(combatmat2),col = colist, col.lim=c(0,1),is.corr = FALSE)
Explore site stability of LICs with ComBat
matchsite <- function(mat, threshold=0.5){
mat[abs(mat)<threshold]<-0
varnames=rownames(mat)
mat = mat[apply(mat, 1, function(s)all(abs(s<threshold)))==FALSE,]
mat = mat[,apply(mat, 2, function(s)all(abs(s<threshold)))==FALSE]
rerowed.mat = rerow(t(mat))
pp=min(nrow(rerowed.mat),ncol(rerowed.mat));
tmpdat=data.frame(col1=varnames) %>%
left_join(.,data.frame(col1=colnames(rerowed.mat)[1:pp],col2=row.names(rerowed.mat)[1:pp], col3=diag(rerowed.mat[1:pp,1:pp])))
tmpdat$col3[ tmpdat$col3==0]<-NA
return(tmpdat)
}
figs = lapply(1:5,function(jj)matchsite((cor_cb[,grep(paste('Site',jj,'.',sep=''), colnames(cor_cb))])))
## [1] "20 row has all zeros, delete this column."
## [1] "17 row has all zeros, delete this column."
## [1] "20 row has all zeros, delete this column."
## [1] "30 row has all zeros, delete this column."
## [1] "15 row has all zeros, delete this column."
## [1] "8 row has all zeros, delete this column."
## [1] "21 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "17 row has all zeros, delete this column."
## [1] "18 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
combat.sitematch = figs[[1]] %>% rename(site1=col2, site1.r=col3) %>%
left_join(.,figs[[2]] %>% rename(site2=col2, site2.r=col3)) %>%
left_join(.,figs[[3]] %>% rename(site3=col2, site3.r=col3))%>%
left_join(.,figs[[4]] %>% rename(site4=col2, site4.r=col3))%>%
left_join(.,figs[[5]] %>% rename(site5=col2, site5.r=col3)) %>%
rename(combatics = col1)
combat.sitematch_cor = combat.sitematch %>%
dplyr::select(combatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
left_join(.,
combat.sitematch %>%
dplyr::select(combatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
reshape2::melt(., id='combatics') %>% group_by(combatics) %>%
summarise(num.matched = 5-sum(is.na(value)), mean.r = mean(value, na.rm=TRUE))
)
combat.sitematch_cor %>%
knitr::kable(., digits=3, caption='ComBat, Site stability') %>%
kableExtra::kable_classic(.,full=F)
ComBat, Site stability
|
combatics
|
site1.r
|
site2.r
|
site3.r
|
site4.r
|
site5.r
|
num.matched
|
mean.r
|
|
Combat,All.IC01
|
0.998
|
0.998
|
0.996
|
0.998
|
0.997
|
5
|
0.997
|
|
Combat,All.IC02
|
0.999
|
0.999
|
0.999
|
0.999
|
0.998
|
5
|
0.999
|
|
Combat,All.IC03
|
0.990
|
0.989
|
0.994
|
0.984
|
0.989
|
5
|
0.989
|
|
Combat,All.IC04
|
0.973
|
0.907
|
0.944
|
0.989
|
0.981
|
5
|
0.959
|
|
Combat,All.IC05
|
0.989
|
0.990
|
0.993
|
0.984
|
0.984
|
5
|
0.988
|
|
Combat,All.IC06
|
0.937
|
0.771
|
0.919
|
0.974
|
0.931
|
5
|
0.907
|
|
Combat,All.IC07
|
0.975
|
0.951
|
0.978
|
0.929
|
0.982
|
5
|
0.963
|
|
Combat,All.IC08
|
0.560
|
0.962
|
0.938
|
NA
|
0.543
|
4
|
0.751
|
|
Combat,All.IC09
|
0.899
|
0.572
|
0.570
|
0.576
|
0.813
|
5
|
0.686
|
|
Combat,All.IC10
|
0.923
|
0.902
|
0.885
|
0.903
|
0.875
|
5
|
0.898
|
|
Combat,All.IC11
|
0.852
|
0.930
|
0.924
|
0.972
|
0.709
|
5
|
0.877
|
|
Combat,All.IC12
|
0.579
|
0.926
|
0.931
|
0.812
|
0.730
|
5
|
0.796
|
|
Combat,All.IC13
|
0.996
|
0.995
|
0.933
|
0.961
|
0.894
|
5
|
0.956
|
|
Combat,All.IC14
|
0.943
|
0.800
|
0.758
|
0.935
|
0.897
|
5
|
0.867
|
|
Combat,All.IC15
|
0.907
|
0.780
|
NA
|
0.839
|
0.841
|
4
|
0.842
|
|
Combat,All.IC16
|
0.810
|
0.661
|
0.972
|
0.783
|
0.723
|
5
|
0.790
|
|
Combat,All.IC17
|
0.595
|
NA
|
NA
|
NA
|
NA
|
1
|
0.595
|
|
Combat,All.IC18
|
NA
|
0.914
|
0.867
|
0.967
|
0.961
|
4
|
0.927
|
|
Combat,All.IC19
|
0.778
|
0.577
|
0.828
|
0.922
|
NA
|
4
|
0.776
|
|
Combat,All.IC20
|
0.906
|
0.528
|
0.705
|
NA
|
0.869
|
4
|
0.752
|
|
Combat,All.IC21
|
NA
|
NA
|
0.829
|
0.842
|
0.587
|
3
|
0.753
|
|
Combat,All.IC22
|
0.912
|
0.927
|
0.926
|
NA
|
0.958
|
4
|
0.931
|
|
Combat,All.IC23
|
0.901
|
0.742
|
0.616
|
0.752
|
0.936
|
5
|
0.789
|
|
Combat,All.IC24
|
0.960
|
0.813
|
0.948
|
0.873
|
0.829
|
5
|
0.885
|
|
Combat,All.IC25
|
NA
|
0.735
|
0.777
|
NA
|
NA
|
2
|
0.756
|
|
Combat,All.IC26
|
NA
|
NA
|
NA
|
0.775
|
NA
|
1
|
0.775
|
|
Combat,All.IC27
|
0.915
|
0.892
|
0.993
|
0.689
|
0.986
|
5
|
0.895
|
|
Combat,All.IC28
|
0.962
|
0.930
|
0.934
|
0.514
|
0.590
|
5
|
0.786
|
|
Combat,All.IC29
|
0.861
|
0.931
|
0.814
|
NA
|
0.585
|
4
|
0.798
|
|
Combat,All.IC30
|
0.582
|
0.661
|
NA
|
0.883
|
0.512
|
4
|
0.659
|
|
Combat,All.IC31
|
NA
|
NA
|
NA
|
NA
|
0.888
|
1
|
0.888
|
|
Combat,All.IC32
|
0.897
|
0.939
|
0.930
|
0.899
|
0.894
|
5
|
0.912
|
|
Combat,All.IC33
|
0.913
|
0.957
|
0.844
|
0.926
|
0.919
|
5
|
0.912
|
|
Combat,All.IC34
|
0.754
|
NA
|
0.942
|
0.746
|
0.931
|
4
|
0.843
|
|
Combat,All.IC35
|
0.766
|
0.507
|
0.648
|
NA
|
NA
|
3
|
0.640
|
|
Combat,All.IC36
|
NA
|
NA
|
0.712
|
NA
|
NA
|
1
|
0.712
|
|
Combat,All.IC37
|
NA
|
NA
|
0.917
|
NA
|
NA
|
1
|
0.917
|
Explore site stability of LICs without ComBat
figs_nocombat = lapply(1:5,function(jj)matchsite((cor_no_cb[,grep(paste('Site',jj,'.',sep=''), colnames(cor_no_cb))])))
## [1] "16 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "18 row has all zeros, delete this column."
## [1] "20 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "17 row has all zeros, delete this column."
## [1] "19 row has all zeros, delete this column."
## [1] "22 row has all zeros, delete this column."
## [1] "22 row has all zeros, delete this column."
nocombat.sitematch = figs_nocombat[[1]] %>% rename(site1=col2, site1.r=col3) %>%
left_join(.,figs_nocombat[[2]] %>% rename(site2=col2, site2.r=col3)) %>%
left_join(.,figs_nocombat[[3]] %>% rename(site3=col2, site3.r=col3))%>%
left_join(.,figs_nocombat[[4]] %>% rename(site4=col2, site4.r=col3))%>%
left_join(.,figs_nocombat[[5]] %>% rename(site5=col2, site5.r=col3)) %>%
rename(nocombatics = col1)
npcombat.sitematch_cor = nocombat.sitematch %>%
dplyr::select(nocombatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
left_join(.,
nocombat.sitematch %>%
dplyr::select(nocombatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
reshape2::melt(., id='nocombatics') %>% group_by(nocombatics) %>%
summarise(num.matched = 5-sum(is.na(value)), mean.r = mean(value, na.rm=TRUE))
)
npcombat.sitematch_cor %>%
knitr::kable(., digits=3, caption='noComBat, Site stability') %>%
kableExtra::kable_classic(.,full=F)
noComBat, Site stability
|
nocombatics
|
site1.r
|
site2.r
|
site3.r
|
site4.r
|
site5.r
|
num.matched
|
mean.r
|
|
noCombat.All.IC01
|
0.998
|
0.997
|
0.992
|
0.995
|
0.994
|
5
|
0.995
|
|
noCombat.All.IC02
|
0.999
|
0.999
|
0.999
|
0.999
|
0.999
|
5
|
0.999
|
|
noCombat.All.IC03
|
0.996
|
0.986
|
0.991
|
0.994
|
0.993
|
5
|
0.992
|
|
noCombat.All.IC04
|
0.986
|
0.903
|
0.937
|
0.926
|
0.984
|
5
|
0.947
|
|
noCombat.All.IC05
|
0.843
|
0.953
|
0.812
|
0.934
|
0.922
|
5
|
0.893
|
|
noCombat.All.IC06
|
0.992
|
0.966
|
0.992
|
0.925
|
0.988
|
5
|
0.972
|
|
noCombat.All.IC07
|
0.944
|
0.940
|
0.798
|
0.804
|
0.763
|
5
|
0.850
|
|
noCombat.All.IC08
|
0.875
|
NA
|
0.870
|
0.638
|
0.913
|
4
|
0.824
|
|
noCombat.All.IC09
|
0.977
|
0.903
|
0.925
|
0.972
|
0.977
|
5
|
0.951
|
|
noCombat.All.IC10
|
0.980
|
0.970
|
0.951
|
0.952
|
0.893
|
5
|
0.949
|
|
noCombat.All.IC11
|
0.901
|
0.934
|
0.912
|
0.893
|
0.726
|
5
|
0.873
|
|
noCombat.All.IC12
|
0.802
|
0.855
|
0.583
|
0.869
|
0.868
|
5
|
0.795
|
|
noCombat.All.IC13
|
0.943
|
0.897
|
0.881
|
0.947
|
0.898
|
5
|
0.913
|
|
noCombat.All.IC14
|
0.961
|
0.938
|
0.869
|
0.834
|
NA
|
4
|
0.900
|
|
noCombat.All.IC15
|
0.977
|
NA
|
0.844
|
0.964
|
0.940
|
4
|
0.932
|
|
noCombat.All.IC16
|
NA
|
NA
|
0.747
|
0.554
|
0.509
|
3
|
0.603
|
|
noCombat.All.IC17
|
NA
|
0.539
|
0.543
|
0.634
|
0.611
|
4
|
0.582
|
|
noCombat.All.IC18
|
0.892
|
0.930
|
0.838
|
NA
|
NA
|
3
|
0.887
|
|
noCombat.All.IC19
|
0.991
|
0.757
|
0.951
|
0.872
|
0.924
|
5
|
0.899
|
|
noCombat.All.IC20
|
0.958
|
0.628
|
0.961
|
0.539
|
0.609
|
5
|
0.739
|
|
noCombat.All.IC21
|
0.707
|
0.908
|
0.746
|
NA
|
NA
|
3
|
0.787
|
|
noCombat.All.IC22
|
0.666
|
0.719
|
0.779
|
0.940
|
0.969
|
5
|
0.815
|
|
noCombat.All.IC23
|
0.626
|
0.754
|
0.835
|
0.843
|
0.775
|
5
|
0.767
|
|
noCombat.All.IC24
|
0.834
|
0.821
|
NA
|
NA
|
0.877
|
3
|
0.844
|
|
noCombat.All.IC25
|
NA
|
0.701
|
NA
|
0.916
|
NA
|
2
|
0.808
|
|
noCombat.All.IC26
|
NA
|
NA
|
0.800
|
0.907
|
NA
|
2
|
0.853
|
|
noCombat.All.IC27
|
NA
|
NA
|
0.634
|
NA
|
0.898
|
2
|
0.766
|
|
noCombat.All.IC28
|
0.965
|
NA
|
NA
|
NA
|
NA
|
1
|
0.965
|
|
noCombat.All.IC29
|
NA
|
NA
|
NA
|
NA
|
NA
|
0
|
NaN
|
|
noCombat.All.IC30
|
0.565
|
NA
|
NA
|
NA
|
NA
|
1
|
0.565
|
|
noCombat.All.IC31
|
NA
|
NA
|
NA
|
NA
|
NA
|
0
|
NaN
|
|
noCombat.All.IC32
|
0.961
|
0.821
|
0.945
|
NA
|
0.789
|
4
|
0.879
|
|
noCombat.All.IC33
|
0.641
|
NA
|
NA
|
0.788
|
NA
|
2
|
0.715
|
|
noCombat.All.IC34
|
NA
|
0.878
|
NA
|
NA
|
NA
|
1
|
0.878
|
|
noCombat.All.IC35
|
NA
|
NA
|
0.862
|
NA
|
NA
|
1
|
0.862
|
aa = nocombat.sitematch %>%
dplyr::select(nocombatics, site1, site2, site3, site4, site5) %>%
left_join(., npcombat.sitematch_cor %>% select(nocombatics, num.matched,mean.r)) %>%
dplyr::filter(num.matched == 5 & mean.r>0.8)
figs = lapply(1:5,
function(j){
obj0 = site$loso[[j]][[1]]
obj=obj0[[3]]
tmp=do.call(cbind,lapply(obj, unlist)) %>%
data.frame(.) %>%
set_names(c('CT','AREA','FAb1000','MDb1000','FC'))
W40 = tmp[order(site$loso[[j]][[1]][[15]], decreasing=TRUE),] %>%
mutate(IC=paste('Site',j,'.IC',stringr::str_pad(1:nrow(tmp),2, pad = "0"), sep=''))
W40[,1:5]=abs(W40[,1:5])/apply(abs(W40[,1:5]), 1,sum)
row.names(W40)<-W40$IC
fig = W40[aa[,j+1],] %>%
mutate(IC=1:nrow(aa)) %>%
melt(., id='IC') %>%
ggplot(.,aes(x=reorder(IC, desc(IC)),y=value,fill=variable))+
geom_bar(stat = "identity",position="fill", color="black", width=0.9) +
ggtitle(paste('Site',j)) +
ylab("Modality Contribution (%)") + coord_flip() + scale_y_continuous(labels=scales::percent) +
theme_minimal() + scale_fill_manual(values=c('salmon','salmon3','aquamarine1',
'deepskyblue1','gold')) +
theme(legend.position='none',
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
return(list(fig,W40[aa[,j+1],]))
})
Morig40 = all$Morig40
W40=data.frame(do.call(cbind,lapply(Morig40[[3]], unlist)))[order(all$Morig40[[15]], decreasing=TRUE),] %>%
set_names(c('CT','AREA','FAB1K','MDB1K','FC'))
W40=W40%>%
mutate(IC=paste('noCombat.All.IC',stringr::str_pad(1:nrow(W40),2, pad = "0"), sep=''))
W40[,1:5]=abs(W40[,1:5])/apply(abs(W40[,1:5]), 1,sum)
row.names(W40)<-W40$IC
fig.all = W40[aa[,1],] %>%
melt(., id='IC') %>%
mutate(IC = gsub('noCombat.All.','',IC)) %>%
ggplot(.,aes(x=reorder(IC, desc(IC)),y=value,fill=variable))+
geom_bar(stat = "identity",position="fill", color="black", width=0.9) +
ggtitle(paste('All Sites')) +
xlab('LICs before ComBat')+
ylab("Modality Contribution (%)") + coord_flip() + scale_y_continuous(labels=scales::percent) +
theme_minimal() + scale_fill_manual(values=c('salmon','salmon3','aquamarine1',
'deepskyblue1','gold')) +
theme(legend.position='none')
grid.arrange(fig.all, figs[[1]][[1]], figs[[2]][[1]], figs[[3]][[1]], figs[[4]][[1]], figs[[5]][[1]],
widths=c(1.5,1,1,1,1,1))

#knitr::kable(rbind(figs[[5]][[2]],W40[aa[,1],]),digits=2)
dims = do.call(rbind,lapply(Morig40[[7]],function(obj)dim(obj[[1]])))[,1]
Explore site stability of LICs after ComBat
figs_combat = lapply(1:5,function(jj)matchsite((cor_cb[,grep(paste('Site',jj,'.',sep=''), colnames(cor_cb))])))
## [1] "20 row has all zeros, delete this column."
## [1] "17 row has all zeros, delete this column."
## [1] "20 row has all zeros, delete this column."
## [1] "30 row has all zeros, delete this column."
## [1] "15 row has all zeros, delete this column."
## [1] "8 row has all zeros, delete this column."
## [1] "21 row has all zeros, delete this column."
## [1] "24 row has all zeros, delete this column."
## [1] "17 row has all zeros, delete this column."
## [1] "18 row has all zeros, delete this column."
## [1] "23 row has all zeros, delete this column."
combat.sitematch = figs_combat[[1]] %>% rename(site1=col2, site1.r=col3) %>%
left_join(.,figs_combat[[2]] %>% rename(site2=col2, site2.r=col3)) %>%
left_join(.,figs_combat[[3]] %>% rename(site3=col2, site3.r=col3))%>%
left_join(.,figs_combat[[4]] %>% rename(site4=col2, site4.r=col3))%>%
left_join(.,figs_combat[[5]] %>% rename(site5=col2, site5.r=col3)) %>%
rename(combatics = col1)
combat.sitematch_cor = combat.sitematch %>%
dplyr::select(combatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
left_join(.,
combat.sitematch %>%
dplyr::select(combatics, site1.r, site2.r, site3.r, site4.r, site5.r ) %>%
reshape2::melt(., id='combatics') %>% group_by(combatics) %>%
summarise(mean.r = mean(value, na.rm=TRUE),num.matched = 5-sum(is.na(value)))
)
combat.sitematch_cor %>%
knitr::kable(., digits=3, caption='ComBat, Site stability') %>%
kableExtra::kable_classic(.,full=F)
ComBat, Site stability
|
combatics
|
site1.r
|
site2.r
|
site3.r
|
site4.r
|
site5.r
|
mean.r
|
num.matched
|
|
Combat,All.IC01
|
0.998
|
0.998
|
0.996
|
0.998
|
0.997
|
0.997
|
5
|
|
Combat,All.IC02
|
0.999
|
0.999
|
0.999
|
0.999
|
0.998
|
0.999
|
5
|
|
Combat,All.IC03
|
0.990
|
0.989
|
0.994
|
0.984
|
0.989
|
0.989
|
5
|
|
Combat,All.IC04
|
0.973
|
0.907
|
0.944
|
0.989
|
0.981
|
0.959
|
5
|
|
Combat,All.IC05
|
0.989
|
0.990
|
0.993
|
0.984
|
0.984
|
0.988
|
5
|
|
Combat,All.IC06
|
0.937
|
0.771
|
0.919
|
0.974
|
0.931
|
0.907
|
5
|
|
Combat,All.IC07
|
0.975
|
0.951
|
0.978
|
0.929
|
0.982
|
0.963
|
5
|
|
Combat,All.IC08
|
0.560
|
0.962
|
0.938
|
NA
|
0.543
|
0.751
|
4
|
|
Combat,All.IC09
|
0.899
|
0.572
|
0.570
|
0.576
|
0.813
|
0.686
|
5
|
|
Combat,All.IC10
|
0.923
|
0.902
|
0.885
|
0.903
|
0.875
|
0.898
|
5
|
|
Combat,All.IC11
|
0.852
|
0.930
|
0.924
|
0.972
|
0.709
|
0.877
|
5
|
|
Combat,All.IC12
|
0.579
|
0.926
|
0.931
|
0.812
|
0.730
|
0.796
|
5
|
|
Combat,All.IC13
|
0.996
|
0.995
|
0.933
|
0.961
|
0.894
|
0.956
|
5
|
|
Combat,All.IC14
|
0.943
|
0.800
|
0.758
|
0.935
|
0.897
|
0.867
|
5
|
|
Combat,All.IC15
|
0.907
|
0.780
|
NA
|
0.839
|
0.841
|
0.842
|
4
|
|
Combat,All.IC16
|
0.810
|
0.661
|
0.972
|
0.783
|
0.723
|
0.790
|
5
|
|
Combat,All.IC17
|
0.595
|
NA
|
NA
|
NA
|
NA
|
0.595
|
1
|
|
Combat,All.IC18
|
NA
|
0.914
|
0.867
|
0.967
|
0.961
|
0.927
|
4
|
|
Combat,All.IC19
|
0.778
|
0.577
|
0.828
|
0.922
|
NA
|
0.776
|
4
|
|
Combat,All.IC20
|
0.906
|
0.528
|
0.705
|
NA
|
0.869
|
0.752
|
4
|
|
Combat,All.IC21
|
NA
|
NA
|
0.829
|
0.842
|
0.587
|
0.753
|
3
|
|
Combat,All.IC22
|
0.912
|
0.927
|
0.926
|
NA
|
0.958
|
0.931
|
4
|
|
Combat,All.IC23
|
0.901
|
0.742
|
0.616
|
0.752
|
0.936
|
0.789
|
5
|
|
Combat,All.IC24
|
0.960
|
0.813
|
0.948
|
0.873
|
0.829
|
0.885
|
5
|
|
Combat,All.IC25
|
NA
|
0.735
|
0.777
|
NA
|
NA
|
0.756
|
2
|
|
Combat,All.IC26
|
NA
|
NA
|
NA
|
0.775
|
NA
|
0.775
|
1
|
|
Combat,All.IC27
|
0.915
|
0.892
|
0.993
|
0.689
|
0.986
|
0.895
|
5
|
|
Combat,All.IC28
|
0.962
|
0.930
|
0.934
|
0.514
|
0.590
|
0.786
|
5
|
|
Combat,All.IC29
|
0.861
|
0.931
|
0.814
|
NA
|
0.585
|
0.798
|
4
|
|
Combat,All.IC30
|
0.582
|
0.661
|
NA
|
0.883
|
0.512
|
0.659
|
4
|
|
Combat,All.IC31
|
NA
|
NA
|
NA
|
NA
|
0.888
|
0.888
|
1
|
|
Combat,All.IC32
|
0.897
|
0.939
|
0.930
|
0.899
|
0.894
|
0.912
|
5
|
|
Combat,All.IC33
|
0.913
|
0.957
|
0.844
|
0.926
|
0.919
|
0.912
|
5
|
|
Combat,All.IC34
|
0.754
|
NA
|
0.942
|
0.746
|
0.931
|
0.843
|
4
|
|
Combat,All.IC35
|
0.766
|
0.507
|
0.648
|
NA
|
NA
|
0.640
|
3
|
|
Combat,All.IC36
|
NA
|
NA
|
0.712
|
NA
|
NA
|
0.712
|
1
|
|
Combat,All.IC37
|
NA
|
NA
|
0.917
|
NA
|
NA
|
0.917
|
1
|
aa = combat.sitematch %>%
dplyr::select(combatics, site1, site2, site3, site4, site5) %>%
left_join(., combat.sitematch_cor %>% select(combatics, num.matched,mean.r)) %>%
dplyr::filter(num.matched==5 & mean.r>0.8)
figs = lapply(1:5,
function(j){
obj0 = site$loso.comba[[j]][[1]]
obj=obj0[[3]]
tmp=do.call(cbind,lapply(obj, unlist)) %>%
data.frame(.) %>%
set_names(c('CT','AREA','FAb1000','MDb1000','FC'))
W40 = tmp[order(site$loso.comba[[j]][[1]][[15]], decreasing=TRUE),] %>%
mutate(IC=paste('Site',j,'.IC',stringr::str_pad(1:nrow(tmp),2, pad = "0"), sep=''))
W40[,1:5]=abs(W40[,1:5])/apply(abs(W40[,1:5]), 1,sum)
row.names(W40)<-W40$IC
fig = W40[aa[,j+1],] %>%
mutate(IC=1:nrow(aa)) %>%
melt(., id='IC') %>%
ggplot(.,aes(x=reorder(IC, desc(IC)),y=value,fill=variable))+
geom_bar(stat = "identity",position="fill", color="black", width=0.9) +
ggtitle(paste('Site',j)) +
ylab("Modality Contribution (%)") + coord_flip() + scale_y_continuous(labels=scales::percent) +
theme_minimal() + scale_fill_manual(values=c('salmon','salmon3','aquamarine1',
'deepskyblue1','gold')) +
theme(legend.position='none',
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
return(list(fig,W40[aa[,j+1],]))
})
Morig40 = all$Morig40.combat
W40=data.frame(do.call(cbind,lapply(Morig40[[3]], unlist)))[order(Morig40[[15]], decreasing=TRUE),] %>%
set_names(c('CT','AREA','FAB1K','MDB1K','FC'))
W40=W40%>%
mutate(IC=paste('Combat,All.IC',stringr::str_pad(1:nrow(W40),2, pad = "0"), sep=''))
W40[,1:5]=abs(W40[,1:5])/apply(abs(W40[,1:5]), 1,sum)
row.names(W40)<-W40$IC
fig.all = W40[aa[,1],] %>%
melt(., id='IC') %>%
mutate(IC = gsub('Combat.All.','',IC)) %>%
ggplot(.,aes(x=reorder(IC, desc(IC)),y=value,fill=variable))+
geom_bar(stat = "identity",position="fill", color="black", width=0.9) +
ggtitle(paste('All Sites')) +
xlab('LICs after ComBat')+
ylab("Modality Contribution (%)") + coord_flip() + scale_y_continuous(labels=scales::percent) +
theme_minimal() + scale_fill_manual(values=c('salmon','salmon3','aquamarine1',
'deepskyblue1','gold')) +
theme(legend.position='none')
grid.arrange(fig.all, figs[[1]][[1]], figs[[2]][[1]], figs[[3]][[1]], figs[[4]][[1]], figs[[5]][[1]],
widths=c(1.5,1,1,1,1,1))
